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ABSTRACT 

The polarization modes of the cosmological microwave background are an invaluable 
source of information for cosmology, and a unique window to probe the energy scale 
of inflation. Extracting such information from microwave surveys requires disentan¬ 
gling between foreground emissions and the cosmological signal, which boils down to 
solving a component separation problem. Component separation techniques have been 
widely studied for the recovery of CMB temperature anisotropies but quite rarely for 
the polarization modes. In this case, most component separation techniques make use 
of second-order statistics to discriminate between the various components. More recent 
methods, which rather emphasize on the sparsity of the components in the wavelet 
domain, have been shown to provide low-foreground, full-sky estimate of the CMB 
temperature anisotropies. Building on sparsity, the present paper introduces a new com¬ 
ponent separation technique dubbed PolGMCA (Polarized Generalized Morphological 
Component Analysis), which refines previous work to specifically tackle the estimation 
of the polarized CMB maps: i) it benefits from a recently introduced sparsity-based 
mechanism to cope with partially correlated components, ii) it builds upon estima¬ 
tor aggregation techniques to further yield a better noise contamination/non-Gaussian 
foreground residual trade-off. The PolGMCA algorithm is evaluated on simulations of 
full-sky polarized microwave sky simulations using the Planck Sky Model (PSM), which 
show that the proposed method achieve a precise recovery of the CMB map in polar¬ 
ization with low noise/foreground contamination residuals. It provides improvements 
with respect to standard methods, especially on the galactic center where estimating 
the CMB is challenging. 

Key words. Cosmology : Cosmic Microwave Background, Methods : Data Analysis, 
Methods : Statistical 
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1. Introduction 


The CMB provides a snapshot of the state of the Universe at the time of recombination. It 
therefore carries invalnable information about the infancy of the Universe and its evolution 
to the current state. After a series of full-sky measurements of the microwave sky (COBE 
( Bennett et al.|[r996 l, WMAP ( Bennett]|2013 1), Planck ( Planck Collaboration et ar]|2013[ ) 
has recently released an accurate estimation of the CMB temperature map. Studying the 
statistical properties of the CMB temperature maps already provides a huge amount of 
information about the primordial state of the Universe. Additionally, the CMB map is a 
radiation with linear polarization that is described by the three polarization modes T, E 
and B ( Hu fc White||199T l; these three modes are observed through the Stokes parameters 
T, Q and U. The information carried by CMB polarization is crucial for: i) improving the 
estimation of the cosmological parameters and break some of their degeneracies, ii) give a 
privileged access to probe the energy scale of inflation through the estimation of the scalar- 
to-tensor ratio r and the tensor spectral index n^. Indeed, inflationary models predict the 
presence of tensor perturbations, which would be traced by the presence of specific CMB 
B-modes. 

From the first full-sky observations of CMB temperature anisotropies by Cosmic Back¬ 


ground Explorer (COBE - see (Bennett et al. 19961) to the Planck 2013 results (Planck 


Collaboration et al.|2013 |, a lot of attention has been paid to the estimation of temperature 
anisotropies by component separation and little to the polarization signal. In contrast to 
temperature, polarization anisotropies are more challenging to measure as their level is 
about 5 orders of magnitude lower than the temperature anisotropies level. 


Polarization anisotropies have first been measured by balloon (Masi et al. 20061 and ground- 


based experiments such as CBI ( Readhe^|2004 1. The WMAP experiments more recently 


provided full-sky observations of the polarization anisotropies (Bennett 20131. The BICEP2 


collaboration initially announced the detection of inflationary B-modes and a measurement 
of the scale-to-tensor ratio r = 0.2, which would be in slight tension with the Planck 2013 
and 2015 results ( [Planck Collaboration 2015b I that puts a stronger constraint r < 0.11 
(95% confidence interval). If this discrepancy could have found various cosmological ex¬ 
planations, it was soon advocated a polarized foregrounds origin (see (Flauger et al.|[2M4l 


and (Mortonson & Seljak 20141). This point was later supported by a study of dust po¬ 


larization made available by the Planck consortium ( Planck Collaboration|[2M4 1 and joint 
BICEP2/Planck investigations ( BICEP2/Keck|2015 1. This controversy highlights the cen¬ 
tral role played by foreground cleaning and component separation to provide an accurate 
estimation of polarized CMB maps. 


Component separation methods have mainly focused on temperature analysis lately (Leach 


et al. 2008 Basak & Delabrouille 2012 Bobin et al. 2013a|b 2014b I. Most component 
separation techniques rely on a linear mixture model. According to this model, each obser- 
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vation of the sky Xi at wavelength Ui is modeled as the linear combination of n components 

Vz — 1, * * * , iff, Xi — ^ ^ 5 (^) 

,n 

where stands for the contribution of component Sj in observation Xi and rii models 
instrumental noise. In this model, the data are assumed to be at the same resolution. This 
model can be more conveniently recast in matrix form: 

X = AS + N. (2) 


where each observation Xi is stored in the z-th row of the observation matrix X, S is the 
source matrix, A is the mixing matrix and N models instrumental noise. From a physical 
viewpoint, the j-th column of A, which we denote by a^, contains the electromagnetic 
spectrum of the source Sj. 


Better known as blind source separation (BSS) in statistics and signal processing (Comon 


& Jutten 20101, component separation aims at estimating the mixing matrix A and the 


sources S from the knowledge of the observations X. This is a classical ill-posed inverse 
problem, which admits an infinite number of solutions. Disentangling physically relevant 
components from erroneous estimates further requires imposing some desired properties on 
the components to be retrieved. In the framework of BSS, component separation methods 
generally enforce some diversity or contrast to distinguish between distinct components. 
To the best of our knowledge, all standard component separation methods for polarized 


CMB observations build upon second-order statistics (Eriksen et al. 2004 Bennett 2013 
Kim et ah]|2009 Aumont fc Macfas-Perez||2007 1. 


Lately, building on different grounds, a novel component separation technique, LGMCA 
(Local-Generalized Morphological Component Analysis - see (Bobin et al. 2013a| |) has 
been introduced. In contrast to standard component separation techniques, the LGMCA 
algorithm relies on the sparse representation of galactic foregrounds, which makes it more 
sensitive to the higher-order statistics of these components. Using GMCA, a joint WMAP- 
Planck high quality full-sky CMB map has been reconstructed. Additionally, this map has 


been shown to be free of detectable thermal SZ (Sunyaev-Zel’Dovich emission (Sunyaev 


& Zeldovich 1970||), and with low foreground contamination (Bobin et al. 2014b|. The 


resulting CMB temperature map has been used to analyze the large-scale anomalies of the 
CMB radiation. Thanks to its low level of foregrounds, masking of the galactic center was 
not necessary for CMB large scale studies ( Rassat et al.|[20T4 Ben-David fc Kovetz||2014 
Aiola et al.||2014 Lanusse et al.|[20T4 |. 
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Contributions 

In this paper, building upon the concept of sparse signal representation, we introduce an 
extension of the LGMCA algorithm to compute an accurate estimate of the CMB po¬ 
larized maps with improved noise/non-Gaussian foreground contamination trade-off. In 
contrast to the processing of CMB temperature data, the estimation of the CMB polar¬ 
ization anisotropies is an even more challenging task since the polarization signal is weak 
and dominated by galactic foregrounds and noise contamination. For that purpose, the 
proposed method, dubbed PolGMCA (Polarized CMC A), builds upon two new key mech¬ 
anisms: i) it makes profit of recent advances in blind source separation (BSS) that allow 
for the efficient separation of sparse and partially correlated sources, ii) it builds upon map 
estimator aggregation, which helps providing a better balance between foreground residu¬ 
als and noise contamination. The PolGMCA algorithm is detailed in Section]^ Numerical 
experiments based on the Planck Sky Model are described in Section which demonstrate 
the ability of the proposed method to achieve good separation performances on the galactic 
center as well as at large scales. 


2. Polarized Sparse Component Separation 

2.1. Component separation for polarized data 

Modeling of polarized data: Describing polarized fields can be done in several ways. They 
are generally measured by the Stokes parameters Q and U, which can further be expanded 
in the Harmonic domain as follows: 


Q±iU = E m; 

£.,m 


( 3 ) 


where ± 2 ^m stands for the spin-2 spherical harmonics basis functions. 

Extending the concept of the wavelet representation for polarization vector fields on the 


sphere has been proposed in (Starck et al. 20091 based on the celebrated "a trous" algo¬ 


rithm ( Starck, J.-L. et ^|2006 |. Building a wavelet transform for Q/U polarization fields 
then consist in building J successive smooth approximations of the Q and U maps via a 
convolution with dyadically rescaled versions of the so-called scaling function h. 


Vj = 0, ••• , J; = hj-kQ 

= hj k U, 


( 4 ) 

( 5 ) 


where J stands for the total number of scales. Following Starck, J.-L. et al. (20061, the 


scaling function is built so as to verify isotropy in the pixel domain. The wavelet coefficients 
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at scale j are then defined as the difference of two consecutive approximations: 


o' 

II 

:> 


Q 

II 
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^j+i 

~ s 

Cj + l 


( 6 ) 

( 7 ) 


In practice, these operations boil down to simple multiplications of the Q and U maps 
spherical harmonics coefficients with the scaling filters {hj}j=o^... _j. 

According to the a trous algorithm, the Q and U maps can be reconstructed by a very 
simple summation: 

.7 

Q = Cj + ^ (8) 

i=i 

U = + (9) 

i=i 


The same polarization field can equivalently be described as E and B polarization modes 
( Zaldarriaga||1998 1: 


E — 

£,m 

B = ^ ^ 

£,m 


^ ^ ty {2^im —2^£m) 



£,m 


( 10 ) 

( 11 ) 


where Yfm stands for the standard 0-spin spherical harmonics functions. The E and B fields 
are therefore usual real scalar fields. The expansions of the E and B fields in the harmonics 
domain is related to the Q and U fields as follows: 


20‘tm + -2a7m 
~ 2 

B .2^im —2^£m 

^ o ■ 


( 12 ) 

(13) 


The E and B fields are therefore derived from the Q and U fields by simply applying twice 
the spin lowering operator to Q + ill and the spin raising operator to Q — iU. Following 


(Starck et al. 20091, extending the isotropic and undecimated wavelet transform for E/B 
polarization fields by constructing formal E and B maps as described in Equation [T^ from 
the Q and U maps and then applying the isotropic undecimated wavelet transform for each 
E and B scalar-valued maps independently. 


Standard approaches for polarized CMB map estimation: Performing component separation 
on a polarization vector field can be performed either in the Q/U coordinate system or 
the E/B field. So far, most standard methods use a description of the components based 
on second order statistics. From that respect, the CMB is more naturally described by 
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its E/B power spectra and cross-spectra in the E and B fields. Therefore, the following 
state-of-the-art methods perform in the E/B fields: 

— Internal Linear Combination - ILC: this component separation technique is well-known 


in the astrophysics community (Eriksen et al. 2004 Bennett 2013 Kim et al. 20091. 


Essentially, it estimates a CMB map, whether in temperature or polarization, with min¬ 
imum variance. The Needlet-ILC ( Basak fc Delabrouille|[2M2 1 performs in the wavelet 
domain; it has been applied to the WMAP data and more recently to the Planck PR2 
data ( Planck Collaboration|2015a l to provide estimates of the CMB polarized maps. In 
statistics, this estimation procedure is well known as the BLUE or Best Linear Unbiased 
Estimator. 

— Spectral matching ICA - SMICA: this component separation method enforces the sta¬ 


tistical independence of the sought after components - for more details, see (Cardoso 


et al. 2002 Delabrouille et al. 20031. SMICA assumes that each component, whether 


the CMB or foreground components, can be modeled as a random Gaussian random 
field with unknown power spectrum. In this setting, the components are estimated by 
enforcing the contrast between the spectrum of the components in the harmonic do¬ 


main. An extension of SMICA to polarized data has been proposed in (Aumont & 
Macias-Perez||2007 1. 


The main limitations of these two component separation techniques are twofold: i) galactic 
and extragalactic foregrounds can hardly be modeled as stationary and homogeneous sig¬ 
nals and ii) they both rely on second-order statistics to disentangle the components, which 
is not well-suited to model galactic foregrounds : they are in nature non-stationary and 
non-Gaussian signals. 

More recently, grounded on the concept of sparsity, we introduced a new component sepa¬ 
ration technique coined CMC A ( Bobin et al.|2013a I. This algorithm is based on the sparse 
modeling of the foregrounds in the wavelet domain. In contrast to the CMB field, the 
galactic foregrounds are better described in the data domain, and therefore in the Q and 
U maps. We will therefore rather perform the CMC A algorithm on the observed Q and U 
maps. 


In early 2015, the Planck consortium released the first estimates of the polarized CMB 


maps that have been computed from the Planck PR2 observations (Planck Collaboration 


2015a I. The proposed 4 maps have been calculated using extensions to polarization of the 


component separation methods used to process CMB temperature maps [Planck Collabo-| 


ration et al. (20131. These methods include: 


— Commander: (Eriksen et al. 20081 is a Bayesian parametric estimation method that 


relies on an explicit parameterized modelling of the sky, which includes CMB, syn¬ 
chrotron, and thermal dust emissions. 
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— SEVEM: ( Fernandez-Cobos et al.|2012 l is a template fitting methods that performs in 
the wavelet domain. For that purpose, templates are derived from differences of input 
Planck observations and used as templates in a foreground cleaning procedure. 

— NILC: A version of the Needlet ILC algorithm has been used to derived polarized 
CMB maps. More precisely, the NILC estimates ILC weights in 13 wavelet bands (see 
( Planck Collaboration||20l'^ I). 

— SMICA: A special version of the SMICA algorithm ( Cardoso et al.|[2002 Delabrouille 
et al.|[2003l has been used to process the Planck PR2 data. In contrast to the original 


SMICA algorithm, spectral parameters are fitted for i < 50. For larger multipoles, it is 
an harmonic ILC. 


Component separation precisely amounts to estimating a mixing matrix A, as described in 
Equation]^ or equivalently unmixing coefficients for ILC-based methods. As in (Aumont & 
|Macias-Perez|2007 Basak fc Delabrouille|2012 |, these parameters can be estimated jointly 
for both the E and B fields, which makes perfect sense for the CMB but not necessarily for 
foreground emissions. Indeed, allowing for different mixing matrices for the Q/U or E/B 
allows for more degrees of freedom to disentangle between complex contaminants, which 
are not perfectly described by the linear mixture model of Subsequently, we propose 
performing the GMCA algorithm independently in the Q and U fields. 


2.2. Sparse component separation with GMCA 


In contrast to standard component separation methods in cosmology, GMCA (Bobin et al. 


2013aI relies on a different separation principle: it makes profit of the naturally sparse dis¬ 


tribution of astrophysical components in the wavelet domain. More precisely, the so-called 
sparsity property means that most of the energy content of the components is concentrated 
in a few coefficients. This property is also shared by the CMB since its power spectrum de¬ 
cays rapidly {i.e. P' at large scale and (A at small scales). This entails that the energy of the 
CMB will be mainly concentrated in a small number of wavelet coefficients. Components of 
distinct physical origins are very likely to exhibit different sparsity patterns. The GMCA al¬ 
gorithm assumes that the components are unlikely to share similar high-amplitude wavelet 
coefficients. Enforcing the sparsity of components in the wavelet domain has indeed been 
shown to be an efficient separation procedure to achieve a clean, low-foreground, estimation 
the CMB temperature anisotropies ( Bobin et al.||201^ 2014b|. 

Let $ denotes a wavelet transform. We assume that each source Si can be sparsely repre¬ 
sented in Sj = ctj^. The data X can be written as 


X = Aa^ -f N , 


(14) 
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where a is &n x T matrix whose rows are aj. The GMCA algorithm estimates a mixing 
matrix A that leads to the sparsest sources S. In practice, this is achieved by solving the 
following optimization problem 


mini|lX-Aa$|j^ + A||a|| 

A,a z 


(15) 


The expression ||a|L stands for pseudo-norm of a, which counts the number of non-zero 


entries of a. The term IjAUp = 


^ ^trace(X^X)^ stands for the Frobenius norm. 


The problem in Equation 15 is solved by using an iterative algorithm that optimizes se¬ 
quentially on the sources S = and the mixing matrix A. According to the model in 
Equation the GMCA algorithm needs to be applied on data that share the same resolu¬ 
tion. For that purpose, the Planck frequency channels are first downgraded to a common 
resolution prior to applying GMCA; estimating a high-resolution CMB map estimate is 
then carried out by performing GMCA on various subsets of the Planck data. We refer the 
interested reader to (Bobin et al. [201^1 for more details about the GMCA algorithm. 


Accounting for the variability of the mixture model : In polarization as well as in temperature, 
most foreground emissions, such as the thermal dust and synchrotron emissions, have 
an electromagnetic spectrum that varies across the sky. From a mathematical point of 
view this makes the mixing matrix A vary across pixels. To account for the heterogeneity 
of the mixture, the GMCA algorithm is applied on sky patches in each wavelet band. 
These different steps gave birth to the LGMCA (Local-GMCA) algorithm, which has been 
described in details in ( Bobin et al.||201^ l. 

The LGMCA algorithm has been implemented and evaluated on simulated Planck data in 
(Bobin et al.|201^|. It has been applied to the Planck PRl data leading to the estimation 


of a low-foreground full-sky map of the CMB anisotropies Bobin et al. (2014b I. 


2.3. Sparse reconstruction of the CMB polarization maps 
Separating non Gaussian and partially correlated components 

Whether in polarization or in temperature, foreground removal is an extremely challenging 
task, particularly : 


In the galactic center: in galactic regions, the foregrounds are, by far, the contam¬ 
inants with the largest emissivity. The rapid variation of the foreground emission laws 
at the vicinity of the galactic center makes component separation strenuous, in partic¬ 


ular for synchrotron and thermal dust emissions. This has been addressed in (Bobin 


et al. 2013a I by localizing the mixture model in the framework of GMCA. However, 


distinct components are also prone to share some similarities in the galactic region: 
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different components have high emissivity in the same regions. Subsequently, distinct 
components will exhibit partial correlations in these areas. This phenomenon makes 
their separation much more straining. 


— At large scales: in this case ~ typically for i < 100 - only few spherical harmonics 
are observed. As a consequence, even theoretically decorrelated components exhibit 
some experimental correlation, customarily dubbed chance-correlations. The underlying 
components can therefore be regarded as partially correlated signals. 


In both cases, the presence of partial correlations between the components will highly 
hamper the performances of component separation methods and weaken their ability to 
efficiently remove galactic foregrounds from the cosmological signal. An extension of the 


GMCA algorithm (Bobin et al. 2014a| has been recently introduced in the signal processing 
community to specifically tackle the separation of partially correlated sources. Its main 
assumption is that partial correlation between components will mainly impact a subset of 
the sparse decomposition coefficients. The rationale of this algorithm relies on a weighting 
scheme that aims at penalizing correlated entries, which are the most detrimental for 
separation. This algorithm first boils down to defining some diagonal weight matrix Q 
( Bobin et al.|2014a l, and secondly to substitute the minimization problem in Equation 15 
with the following one: 


min - Trace {(X - Aa$)Q(X - Aa$)^} -b A |ja|L 
A, a 2 


(16) 


In (Bobin et al. 2014aI, it has been shown that partially correlated samples are related 


to non-sparse columns of the source matrix S. Based on this relationship, it adaptively 
updates the weight matrix with respect to the estimated sources during the iterations of 
the algorithm. 

This algorithm has been shown to be very effective to disentangle partially correlated 
sources. For more details, we refer the interested reader to (Bobin et al.||2014a|. 


Improving the noise/foreground residual trade-off by combining estimators 
Whether in temperature or in polarization, the accurate estimate of CMB anisotropies 
requires solving somehow antagonistic signal processing issues: 


— Low noise contamination : polarized CMB data are highly dominated by instrumen¬ 
tal noise, by several orders of magnitude. An effective CMB recovery method should 
therefore be able to deliver polarized CMB map estimates with low noise contamination. 
Since instrumental noise in Planck-like data is generally approximated by a stationary 
Caussian random field. Hence methods based on second-order statistics (SOS) in the 
harmonic domain, such Harmonic ILC (HILC - ( Kim et aL]|2009[ |), are very efficient at 
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small scales where the noise is dominating. 


— Low foreground residuals : foregrounds, and especially galactic foregrounds, are non- 
Gaussian and non-stationary signals. As a consequence, component separation methods 
that rely on SOS, and particularly in the Harmonic domain, are not well suited to 
accurately remove foreground components. In ( Bobin et al.j2013a l, we emphasized that 
the use of sparsity in spatially localized signal representations, such as wavelets, provides 
an effective strategy to estimate a low-foreground CMB map. 

These two points underlines that component separation for CMB polarization data differs 
from the case of temperature where the noise is not as large compared to the level of 
foreground components. Neither SOS-based nor sparsity-based techniques are perfect can¬ 
didates to estimate low-foreground and low-noise CMB polarization maps. Subsequently, 
we rather propose to combine estimates delivered by SOS-based and sparsity-based com¬ 
ponent separation methods. This aggregation of complementary estimators is inspired by 
advances in statistics and estimation theory ( Yang|[20M l. 

For the sake of generality, we assume that we have access to P unbiased estimators of the 
CMB {ci}i=i^... ^p, which basically differ from each other by their noise contamination and 
level of foreground residuals. The goal of estimator aggregation is to find some weights 
- ,p so that the combined estimator is expressed as a linear combination of these 
estimators: 


c = ^ 7r,Ci 

Z=1 


(17) 


The resulting estimator is guaranteed to be unbiased as long as the weights sum up to unity: 

Furthermore, the weights are estimated by minimizing some cost function 

77.(7ri, • • • ,7rp) : 


min 77.(7ri,-'' ,7rp) s.t. 

■ ,p 


= 1 . 


(18) 


which enforces some desired properties for the combined estimator. 

Essentially, the proposed estimator aggregation procedure will amount to combining a SOS- 
based CMB map estimates, delivered by the HILC, with sparsity-based estimates provided 
by the LGMCA algorithm. The LGMCA algorithm is well suited for accurately removing 
foreground residuals but does not provide a low-noise estimate. To that respect, estimator 
combination is mainly implemented so as to minimize the impact of instrumental noise. 
Since noise is best modeled in the Harmonic domain with second-order statistics, it is 
therefore natural to combine LGMCA-based estimates with HILC by seeking the weights 
- ,p that minimize the variance of the estimated CMB in the harmonic domain. 
In practice, deriving a combined estimator in the harmonic domain is done by decomposing 
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the harmonic domain in overlapping bins of multipoles of fixed width A, which will be 
indexed by 1. A set of weights {TTij}i-i^... is then estimated in each band of multipoles. 
Consequently, the aggregation of the P estimators is carried out by solving the 

following minimization problem in each band I of multipoles: 


mm 

}i=l, - ,p 




S.t. 




(19) 


i2 


i=l 


where |j . denotes the Euclidean norm, the term Ci^i stands for the ^-th harmonic multi¬ 
pole of the i-th CMB map, and the values bi i are the the beam of different maps at the l-th 
harmonic multipole. This estimator is better known as the Best Linear Unbiased Estimator 
(BLUE) in the statistics community or the ILC in astrophysics. The different weights are 
then re-combined to provide values for each multipole £. 

Interestingly, the LGMCA algorithm amounts to applying the GMCA algorithm on vari¬ 
ous subsets of observations. This yields various estimates of the CMB maps with various 
spatial resolutions : estimated maps with the lowest resolution are computed from a larger 


number of observations and vice-versa. We have pointed out in (Bobin et al. 2014bI that 


applying component separation to a large number of observations allows for more degrees 
of freedom to better provide a clean estimate of the CMB map. The current version of 
the LGMCA algorithm performs by using a single CMB map estimate in each wavelet 
band. From the perspective of estimator aggregation, the LGMCA algorithm already com¬ 
bines complementary CMB maps but in a rather suboptimal way. Again, the framework 
of estimator aggregation can help providing a more efficient combination of GMCA-based 
estimators with a better balance between noise contamination (he. using maps with dif¬ 
ferent spatial resolutions) and foreground residual (i.e. using maps obtained from various 
subsets of data). 

The PolGMCA algorithm therefore i) derives 5 estimates of the CMB map from various 
subsets of observations and with different spatial resolutions (see Section |3.1| ), ii) com¬ 
bines these sparsity-based estimates together with the HILC estimates, and iii) it provides 
an accurate estimation of large-scale multipoles without being highly sensitive to chance- 
correlations. Conversely, since the estimator aggregation procedure is based on second-order 
statistics, it may be prone to chance-correlations at large-scale. The aggregation weights 
are therefore set to zero for £ < 100 with the exception of the LGMCA low-resolution 
estimator. 

To illustrate how estimator aggregation performs. Figure displays the values of the ag¬ 
gregation weights for each of the 5 Q map estimates computed with LGMCA and the 
HILC map. These results have been obtained by performing the PolGMCA algorithm to 
the Planck-like PSM simulations described in Section [S] 

It is very important to notice that, as expected, the contribution of each of the maps de- 
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livered by the LGMCA algorithm depend directly on their resolution: low-resolution maps 
have significant weights at large-scales and vice-versa. Moreover, HILC mainly contributes 
to the final estimate at small scales {£ > 700) and is predominant at the smallest scales 
(£ > 1500) where noise is very likely the predominant source of contamination. The global 
evolution of these weights clearly corroborates the original motivation of estimator aggre¬ 
gation : i) give more weight to sparsity-based estimators when non-Gaussian foregrounds 
are likely the most prominent, and ii) give more weight to the SOS-based estimator when 
instrumental noise is the main contamination. 


Weights estimated for estimator aggregation 



Fig. 1. Weights computed by estimator aggregation. 


3. Numerical experiments 

3.1. Simulations and experimental protocol 


Simulations By the end of 2015, the Planck satellite will provide the finest observations of 
the polarized sky in the frequency range 30 - 353GHz. We naturally propose evaluating the 
proposed component separation method with simulated polarized Planck data. For that 
purpose, we will make use of simulations generated with the Planck Sky Model (PSM)[^in 


(Delabrouille et al. (20131). 


PSM simulations The PSM includes the best of the current knowledge of most polarized 
astrophysical signals and foregrounds, as well as simulated instrumental noise and beams. 
In detail, the simulations were obtained using the publicly available PSM version 1.7.8 with 
the following: 

— Frequency channels: the simulated data are composed of the 7 LFI and HFI chan¬ 
nels at frequency 33,44, 70,100, 143, 217, 353 GHz. The frequency-dependent beams are 
perfect isotropic Gaussian PSF with FWHM ranging from 5 arcmin at 217, 353 GHz to 
33 arcmin at 33GHz. These simulations do not account for the effect of the scanning 
strategy adopted by full-sky surveys; indeed, these survey strategies are akin to produce 
non-isotropic and spatially varying beams, which are more elongated in the scanning 

^ For more details about the PSM, we invite the reader to visit the PSM website: http://www. 
ape.univ-paris7.fr/~delabrou/PSM/psm.html 
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directions. 

— Instrumental noise: the sky coverage is not uniform and homogeneous, since some 
areas are more often observed than others. The noise statistics (i.e. noise variance map 
for each frequency channel) are assumed to be known accurately and consistent with 
the Planck noise variance maps. The simulations account for the inhomogeneity of the 
noise statistics but not for its correlation along the scanning direction. Noise will be 
assumed to be uncorrelated. 

— Cosmic microwave background: the CMB map is a correlated random Gaussian re¬ 
alization entirely characterized by its T, E, B power spectrum and T-E cross-spectrum. 
In the simulations, the CMB intensity and polarization power spectra and cross-spectra 
are generated according to the A-CDM model derived from the Planck 2013 results. The 
simulated CMB is Gaussian, and no non-Gaussianity (e.g. leasing, ISW, fAri) has been 
added. The simulated polarized CMB maps are free of B-mode. 

— Point sources: infrared and radio sources were added based on existing catalogs at 
that time. In the PSM, polarized point sources are attributed a random polarization 
degree and uniformly random polarization angle. 


The galactic polarization model, which are composed of the dust and synchrotron emis¬ 
sions. These emissions are based on the PSM mamd2008 model IMiville-Deschenes et al.l 


(20081 that is parameterized by the intrinsic polarisation fraction. By default, it is set to 


15% to match the WMAP data. 
More precisely. 


— Polarized dust emissions: The polarization of galactic thermal dust emission is due 
to the partial alignment of elongated dust grains with the Galactic magnetic field. It is 
however poorly known; the current PSM simulations are extrapolation from the thermal 
dust intensity. The emission of the thermal dust is modeled by a gray-body emission 
law. PSM simulations has been shown to be in good agreement with the Archeops ob¬ 
servations of polarized dust (Ponthieu et al.|[2005| and WMAP 94GHz measurement. 


Polarized synchrotron emission: The polarized synchrotron emission is obtained 
by extrapolating WMAP 23GHz Q and U maps using the same power law model as 
for the temperature map. This approach has shown to be consistent with the WMAP 
polarized synchrotron maps (Gold et al.||2011 1. 
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Band 

^ common resolution 

# sources 

Nominal patch size 

1 

33 arcmin 

7 

N/A 

2 

24 arcmin 

6 

64 

3 

14 arcmin 

5 

64 

4 

10 arcmin 

4 

32 

5 

5 arcmin 

3 

16 


Table 1. Parameters of the local and multiscale mixture model used in LGMCA and PolGMCA. 


Comparison protocol So far, the component separation methods used to estimate CMB 
polarization maps are limited to ILC-based techniques (Basak fc Delabrouille|20f2| and an 


extension of SMICA to polarization (Aumont & Macias-Perez 20071. Contrary to LGMCA, 


none of these codes is publicly available. The recent analysis of the Planck PRl tempera¬ 


ture data (Planck Collaboration et al. 20131 revealed that ILC-based component separation 


techniques, and especially its harmonic variant, provided some of the best results in the 
noise-dominant regime: the SMICA map has been obtained by applying Harmonic ILC 
(HILC) for ^ > 1500. Hence, it is very likely that similar results will hold true in the case 
of polarization. Consequently, this section will focus on comparing the CMB map derived 
from HILC (described in the appendix]^ and the new PolGMCA method, which combines 
HILC with the updated version of the LGMCA algorithm. 


The LGMCA algorithm requires defining four parameters: i) the number of sources is 
set to be equal to the number of channels, ii) the common resolution, and iii) the nominal 


patch size. All these parameters are given in Table 3.1 and are chosen as in Bobin et al. 


( |2014b| . 


3.2. CMB map estimates 

Figures and display the estimated CMB Q and U maps, which have been computed 
with HILC and PolGMCA. The maps have been smoothed at the resolution of 1 degree so 
as to better highlight the large-scale residuals features and remove high-frequency noise, 
thus easing visual inspection. These figures first show that the proposed PolGMCA al¬ 
gorithm achieves a much better cleaning of the galactic center; the HILC map clearly 
contain significant large-scale foreground residual, which are mainly limited to the galactic 
center. However, while showing less residuals on the galactic center, the PolGMCA map 
seems to be slightly more contaminated by instrumental noise as seen by seemingly larger 
background fluctuations. To a lower extent, the same conclusions can be drawn from the 
estimated U maps features on Figure]^ 

Figure 1^ shows the total foreground residuals which contaminates the estimated CMB Q 
maps. These maps are obtained by applying the unmixing parameters to the input fore¬ 
grounds. These maps have further been smoothed to a common resolution of 1 degree. Both 
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HILC and PolGMCA visually yield low levels of foreground residuals at large galactic lati¬ 
tudes. As noticed earlier in this paragraph, the PolGMCA algorithm achieves a significantly 
better foreground cleaning at the vicinity of the galactic center. 


Input CMB Q map 


CMB Q map - PolGMCA 


-5.0 


5 0 fiK 


■ ■ ■ 

-5.0 - 




5.0 /iK 





Fig. 2. Estimated CMB Q maps at 1 degree resolution: Input CMB (left), HILC (middle), 
PolGMCA (right). 


Input CMB U map 


CMB U map HILC 


CMB U map PolGMCA 



Fig. 3. Estimated CMB U maps at 1 degree resolution: Input CMB (left), HILC (middle), 
PolGMCA (right). 



Total residual Q map - HILC 


Total residual Q map — PolGMCA 







Fig. 4. Foreground residual in the estimated Q maps at 1 degree resolution: Input CMB 
(left), HILC (middle), PolGMCA (right). 


Input CMB U map 


Total residual U map HILC 


Total residual U map PolGMCA 








Fig. 5. Foreground residual in the estimated U maps at 1 degree resolution: Input 
CMB (left), HILC (middle), PolGMCA (right). 


Total foreground residual In this section, we evaluate the quality of estimation of the polar¬ 
ized CMB E maps, which we computed with HILC and PolGMCA. It goes without saying 
that the purpose of the article is to provide a component separation methods that provides 
clean, low-foreground, polarized CMB maps at high galactic latitudes as well as on the 
galactic center. For that purpose, the power spectra of the individual and total foreground 
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residuals have been computed from various values of the sky coverage -25%, 55% and 85% 
- and without any masking (full-sky). These masks are common to all the methods; they 
have been computed based on the foreground residuals of both the two maps, so as not to 
favor/penalize one of these maps. Since these comparisons are based on simulations, the 
noise and CMB angular power spectra are the actual spectra of the input simulation; they 
are therefore not based on Monte-Carlo simulations. Furthermore, these angular power 
spectra are display with Scr error bars based on 100 noise simulations. The instrumental 
noise being one of the major contaminant of the estimated polarized CMB maps, the un¬ 
certainties with respect to noise allow to assess quite accurately the statistical significance 
of the foregrounds residuals. 

Figure displays the power spectrum of the total foreground residual (i.e. synchrotron, 
thermal dust and point sources) for HILC and PolGMCA for various values of the sky 
coverage. At large-scale, for I < 100, the PolGMCA algorithm provides a lower level of 
foreground residuals whether it is at large galactic latitudes (sky coverage of 25%) or on 
the galactic center (full-sky estimate). This enhancement is very likely the consequence of 
the new separation mechanism described in Section |2.3| which makes the estimation of the 
E CMB map less sensitive of chance-correlations, but at the cost of an increase of the noise 
level. 

At intermediate scales - for 100 < i < 600 - PolGMCA and HILC provide very similar 
results on 55% of the sky. HILC seem to perform slightly better for high galactic latitudes. 
Both methods yield estimates with very low foreground residuals at high galactic latitudes, 
which are about one order of magnitude lower than both the noise and the cosmological 
signal. For larger sky coverage, the level of foreground residuals rapidly increases in the 
CMB E map delivered by the HILC algorithm while PolGMCA produces a clean full-sky 
estimate of the E map at intermediate scale. At smaller scales - i > 600 - PolGMCA 
performs better than HILC with similar noise level but lower foreground residuals. Inter¬ 
estingly, at very large-scales - £ < 10 - the level of foregrounds is significantly lower than 
the level of CMB; this suggests that large-scale anomalies in polarization could also be 
carried out without masking, similarly to what has been done for the temperature large 
scale anomalies studies (Rassat et al. ]M4| . This result makes PolGMCA a good candidate 
to analyze the large-scale anomalies of the CMB map. 


Figure [^displays 3cr error bars based on 100 noise simulations, which first reveal that 
outside the galactic (i.e. for masks larger than 85%) the discrepancy between these two 
methods is not significant with respect to noise uncertainty. In contrast, when no mask 
is applied, the improvements led by the PolGMCA algorithm are statistically significant. 
This highlights the very good performances of the PolGMCA algorithm at the vicinity of 
the galactic center. 
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TOTAL FRG power spectrum - E map - MASK Isky 25 

- Harmonic ILC 

PoIGMCA 
- InputCMB 



TOTAL FRG power spectrum - E map - MASK Isky 85 


- Harmonic ILC 

PoIGMCA 



TOTAL FRG power spectrum - E map - MASK Isky 55 


- Harmonic ILC 

PoIGMCA 
- InputCMB 



TOTAL FRG power spectrum - E map - NOMASK 


- Harmonic ILC 



Fig. 6. Total foregrounds - top-left : power spectrum of total foregrounds residual for the E 
map from a sky coverage of 25%. top-right : from a sky coverage of 55%. bottom-left : from a 
sky coverage of 85%. bottom-right : full-sky. 


Individual foreground residuals This paragraph investigates the contamination of the polar¬ 
ized CMB map estimates, and especially the CMB E map from individual foregrounds. 


— Thermal dust Figure [^features the power spectrum of the residual thermal dust con¬ 
tamination for various sky coverage ranging from 25% to full sky. For intermediate and 
small scales {i.e. £ > 1000), thermal dust is by far the major contaminant for both the 
HILC and PoIGMCA E maps, whether it is on the galactic center or at high galactic 
latitudes. In this range of multipoles, HILC performs slightly better than PoIGMCA 
for sky coverage smaller than 55%. For larger sky coverage, the ability of PoIGMCA to 
deliver a clean estimate of the CMB map on the galactic center greatly helps limiting 
the increase of thermal dust residual from areas close to the galactic center. It is im¬ 
portant to notice that the E map delivered by the PoIGMCA algorithm shows a level 
of thermal dust residual that is about one order of magnitude lower than the level of 
CMB for £ < 1200, even for a sky coverage as large as 85%. 

At small scales, for £ < 100, the CMB E map computed with PoIGMCA contains sig¬ 
nificantly lower thermal dust residuals at all latitudes. HILC provides a slightly lower 
dust level in the range 100 < £ < 1000 at high galactic latitudes. 

— Synchrotron the power spectrum of the synchrotron residual in the HILC and PoIGMCA 
E maps are displayed in Figure HILC a provides lower synchrotron contamination 
for £ < 20, whether it is on the galactic center or at larger galactic latitudes. However, 
PoIGMCA provides significantly lower synchrotron residuals for £ > 100, even for high 
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galactic latitudes. 

— Point sources Figureshows the power spectrum of the point sources contamination. At 
large-scale, PolGMCA exhibits a slightly larger amount of point sources residual, for all 
sky coverage. PolGMCA performs slightly better than HILC at intermediate and small 
scales. It is interesting to point out that that for large sky coverage, the E map obtained 
with PolGMCA shows a significantly lower level of point sources residual for I > 100; 
this highlights the ability of the PolGMCA algorithm to remove highly non-Gaussian 
and non-stationary components like point sources, especially on the galactic center. For 
astrophysical studies, the most prominent point sources will be masked to limit their 
impact. The level of point sources residual is therefore expected to be significantly lower 
for both component separation techniques in this setting. 


IrermaldusI power spectrum - 




- Harmonic 
PolGMCA 

- InputCMB 



thermaldust power spectrum - E map - MASK fsky 85 


- Harmonic 
PolGMCA 



- Harmonic 
PolGMCA 

- InputCMB 



thermaldust power spectrum - E map - NOMASK 



Fig. 7. Thermal dust - top-left : power spectrum of the thermal dust residual in the estimated 
E map from a sky coverage of 25%. top-right : from a sky coverage of 55%. bottom-left : from 
a sky coverage of 85%. bottom-right : full-sky. 


3.3. Towards an accurate full-sky estimate of the CMB polarization maps 


In this section, we investigate with more precision the quality of estimation of the GMB 
polarization map across the sky. Indeed, one of the objectives of this article is to allow the 
accurate full-sky estimation of polarized GMB maps. For that purpose, and following an 
analysis method we have used (Bobin et al. |2013a|b| |2014b| |, the accuracy of GMB map can 
be evaluated by computing the level of foreground residuals in bands of latitudes, in the 
wavelet domain. Using wavelets, one can analyze the estimated GMB maps with a good 
localization in space as well as in harmonic space. 
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synchrotron power spectrum - E map - MASK fsky 25 


- Harmonic ILC 

PoIGMCA 
- InputCMB 



synchrotron power spectrum - E map - MASK fsky 85 


- Harmonic ILC 

PoIGMCA 



synchrotron power spectrum - E map - MASK fsky 55 


- Harmonic ILC 

PoIGMCA 
- InputCMB 


synchrotron power spectrum - E map - NOMASK 



- Harmonic ILC 

PoIGMCA 



Fig. 8. Synchrotron - top-left : power spectrum of the synchrotron residual in the estimated 
E map from a sky coverage of 25%. top-right : from a sky coverage of 55%. bottom-left : from 
a sky coverage of 85%. bottom-right : full-sky. 


ps power spectrum - E map - MASK fsky 25 

- Hormonlc ILC 

PoIGMCA 
- InputCMB 



ps power spectrum - E map - MASK fsky 85 



ps power spectrum - E map - MASK fsky 55 

- Hormonlc ILC 

PoIGMCA 
- InputCMB 



s power spectrum - E map - NOMASK 



Fig. 9. Point sources - top-left : power spectrum of the point sources residual in the estimated 
E map from a sky coverage of 25%. top-right : from a sky coverage of 55%. bottom-left : from 
a sky coverage of 85%. bottom-right : full-sky. 


More precisely, the estimated Q and U maps have been decomposed in 6 wavelet bands 


see (Starck, J.-L. et al. 20061. Each of the 6 bands is then split uniformly into 18 bands 
of latitudes; each band of latitude has width 10°. Figure [T0| reports the level of energy 
(Euclidean norm) of the total foreground residuals (i.e. thermal dust, synchrotron, and 
point sources) for the 6 wavelet scales. By convention, the wavelet scale 0 corresponds to 
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the smallest scales > 1600) and scale 5 to the large scales (50 < i < 100). 

It is first interesting to notice that, at all scales, HILC has a level of foreground residuals 
that exceed the level of CMB on the galactic center, in latitudes of about ±15°. This is 
precisely at the vicinity of the galactic center that the stationarity assumptions, made in 
the framework of HILC, is probably the least valid and where the use of spatial localiza¬ 
tion and sparsity is the most beneficial for component separation. Consequently, the CMB 
map estimated by the PolGMCA algorithm has significantly lower foreground residuals 
for galactic latitudes lower than 15°. More precisely, in the galactic plane, HILC Q map 
has an order of magnitude higher foreground residuals than PolGMCA, in wavelet bands 
corresponding to ^ > 400. 

For the first three wavelet scales, the PolGMCA Q and U maps have a level of foreground 
residuals which is constantly lower than the level of CMB all over the sky. At larger scales, 
CMB is dominated by the foregrounds only for latitudes which are lower than ±10°. 

For galactic latitudes larger than ±25°, these plots confirm the very good efficiency of 
the HILC algorithm in regions where instrumental noise is, by far, the most prominent 
contaminant in the data. The PolGMCA algorithm generally provide similar or slightly 
higher levels of foregrounds in high galactic regions, but still much lower than the level of 
instrumental noise. 

These numerical results highlight the efficiency of the proposed PolGMCA algorithm to 
provide a very accurate estimation of the polarized CMB maps at the vicinity of the galac¬ 
tic center while providing residual levels that are competitive to the very effective HILC 
algorithm at high galactic latitudes at intermediate and small scales. The PolGMCA is 
shown to be particularly effective at large-scales {i.e. for £ < 200), where it yields esti¬ 
mates of the CMB Q and U maps with significantly lower levels of total foreground at all 
galactic latitudes. 


4. Conclusion 

We present a novel component separation method which is dedicated to the estimation of 
the CMB polarization modes. This approach builds upon the GMCA algorithm, which em¬ 
phasizes on the sparsity of the sought after sources in the wavelet domain to separate them 
out. The proposed dubbed PolGMCA extends the basic GMCA algorithm with two signif¬ 
icant improvements: i) it allows for the efficient separation of partially correlated sources 
thanks to recent improvements of the GMCA algorithm, and ii) it aggregates sparse and 
second-order based CMB estimators to benefit from their complementary characteristics 
in terms of propagated noise and foreground residuals. Subsequently, we show that the 
proposed PolGMCA algorithm exhibit the following interesting properties : 
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Total loreground MSE in wavelet band MO - Q map - 1536 < I < 3072 


Total foreground MSE rn wavelet band #1 - Q map - 768 <l < 1536 





Fig. 10. Total foregrounds residuals in the Stokes parameter Q map: mean squared 
error (MSE) per band of latitude in the wavelet domain. Dashed lines represent the noise level for 
each component separation method. 


— It provides an improved estimation of the polarized CMB maps particnlarly at large- 
scales, typically for i < 200, outside the galactic center, where the ability to account for 
the partial correlation the components is the essential. However, these improvements 
are not significant with respect to the noise-related uncertainty. 

— It provides very accurate estimates of the polarized CMB maps in regions where fore¬ 
grounds, which are by nature non-Gaussian and non-stationary, dominate. This is specif¬ 
ically the case at the vicinity of the galactic center at latitudes ±15°. In this setting, 
the PolGMCA algorithm benefits to a large amount of its ability to account for the 
partial correlation of the components to be removed 

— In noise-dominated regions of the sky, the polarized CMB maps provided by the Pol¬ 
GMCA algorithm have levels of foregrounds which are close to the HILC algorithm, 
which is assumed to be close to optimality in this noise regime. 


Altogether, these results make the PolGMCA algorithm a good candidate for galactic 
studies in polarization, studying large-scale anomalies and testing the detection of B-modes 
from Planck data. 
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Total foreground MSE in wavelet Band #0 - U map - 1536 < I < 3072 


Total foreground MSE in wavelet band - U map - 768 < I < 1536 




Harmonic ILC - 

PoiGMCA 


InputCMB - 


Harmonic ILC - 

PoiGMCA 


InputCMB - 


Total foreground MSE in wavelet bi 


Total foreground MSE in wavelet band #5 - U map - 48 <1 < 96 




Fig. 11. Total foregrounds residuals in the Stokes parameter U map: mean squared 
error (MSE) per band of latitude in the wavelet domain. Dashed lines represent the noise level for 
each component separation method. 
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Appendix A: Harmonic ILC implementation 

Harmonic ILC has been widely used, mainly to estimate CMB maps from WMAP data in 
temperature data (Tegmark et al.|2003l and polarization (Kim et al.|2009 1. More recently, 


the SMICA CMB polarized maps derived from the Planck PR2 data (Planck Collaboration 


2015a I are equivalent to a HILC solution at high i multipoles. The appeal of HILC resides 


in its ability to precisely account for the different beams at each channel as well as its 
accuracy to model stationary Gaussian random fields like the CMB. 

In this paper, and inspired by ( Tegmark et al.|2003 Kim et al.|2009 1, we have implemented 
a HILC method that performs into two steps: 
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— Fitting the ILC weights: ILC estimates weights to obtain a minimum variance solution. 
In the case of HILC, this step is performed independently on overlapping bands of 
multipoles that span the desired range of multipoles [0,2500]. The bands have a fixed 
width A = 80; this width has been tuned so as to provide the best estimation accuracy 
in PSM simulations, as a trade-off between good foreground modeling (narrow bands) 
and statistical robustness (broad bands). All the 81 bands used to derive ILC weights 
are displayed on the left panel of Figure jATT] Let be the beam of the f-th frequency 
channel, be the filter defining the Lth band, and be the coefficient of 

the z-th frequency channel, the weights ,m in the l-th band are computed by 

minimizing the following problem: 



Fig. A.l. Bands of multipoles used to derive HILC weights: left All the 81 overlapping 
bands used to compute the ILC weights in the harmonic domain. Right: examples of a 5 different 
bands for better visualization. 


— Deriving (.-dependent weights: from the 81 bands, one derives weight per bands and 
observations. A simple interpolation procedure is used so as to derive a weight per ( 
and per observation i. This interpolation step is performed to provide weight vectors 
TT* = [7ri(0), • • • ,Tr,{(),--- ,7ri(fmax)] as follows: 


Vz = 1, 




Ei wi,e 


(A.2) 


which is essentially a simple interpolation by weighted averaging, where the filters wi 
are used as weights. 


The exact same HILC algorithm is used to aggregate the CMB estimators in the PolGMCA 
algorithm described in 
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